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So.  Introduction  and  Summary 

A  variety  of  Markov  chains  in  continuous  tine  modeling  stochastic 
systems  of  applied  interest  have  for  their  state  space  a  rectangular  lat¬ 
tice  of  states  8  =  {(j,n)  :  0  j  <  J ,  O^njtN}.  When  the  number  of 
states  (J+1)(N+1)  is  large,  say  >  100,  evaluation  numerically  of  the 
ergodic  distribution,  and  moments  of  exit  times  and  entrance  times  to  sub¬ 
sets  of  interest  is  costly  and  simulation  is  often  resorted  to. 

For  many  such  chains,  changes  in  column  index  j  or  row  index  n  at 
transition  epochs  have  values  0,  ♦  1.  The  chains  may  then  be  described 
as  column- continuous  and  row-continuous  respectively. 

When  such  row- continuity  is  present,  for  example, 

systematic  treatment  of  the  row  subsets  of  states  as  probabilistic  entities 
provides  a  theoretical  basis  for  the  discussion  of  the  chain,  and  algo¬ 
rithms  for  the  description  of  the  chain  involving  matrices  of  order  J+l 
rather  than  (J+1) (N+l) ,  better  suited  to  the  capacity  constraints  of  com¬ 
puters.  The  procedure  may  therefore  be  described  as  rank  reducing. 

Algorithms  based  on  such  treatments  of  rows  as  entities  have  been 
developed  by  M.  Neuts  [10, 11],  when  N  is  infinite,  for  the  study  of  queues 
with  service  times  or  interarTival  times  describable  in  terms  of  "phases". 
His  algorithms  deal  with  state  spaces  for  which  N*«»  and  the  transition 
rates  for  the  chain  are  independent  of  row  index  n,  except  near  the  boundary 
n*0.  His  methods  are  oriented  largely  toward  the  ergodic  behavior  of  such 
chains . 

The  present  treatment  is  primarily  directed  towards  finite  markov 
chains  with  transition  rates  dependent  on  both  j  and  n.  Entrance  and 
exit  time  moments  are  obtained,  along  with  the  ergodic  distributions. 


Entry  and  exit  time  distributions,  obtained  via  the  Laguerre  transform 
[8],  will  be  discussed  elsewhere. 

It  should  be  emphasized  that  the  row  or  column  orientation,  natural 
for  some  systems,  may  be  an  effective  tool  for  chain  description  even 
when  no  natural  row  or  column  meaning  is  present. 

In  the  first  section,  the  basic  bivariate  process  is  described  and 
notation  is  developed.  Several  motivating  examples  are  given.  Subsequent 
sections  develop  the  methodology,  and  algorithmic  procedures,  and  discuss 
computer  efficiency.  In  a  concluding  section  a  tandem  queue  with  Poisson 
arrivals,  exponential  service  of  different  rates,  multiple  servers,  and 
finite  waiting  rooms  is  presented. 
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§1.  The  Bivariate  Markov  Process  B(t)  »  [J(t),N(t)] 

Consider  a  bivariate  Markov  process  B(t)  *  [J(t),N(t)]  on  8  »  JxA/ 
where  J  «  {j:  0  <  j  _<  J},  W  *  {n:  0  <  n  <  N}.  In  a  typical  context,  the 
process  J(t)  is  a  finite  Markov  chain  in  continuous  time  (independent  of 
N(t)),  but  N(t)  is  not  Markov  and  depends  on  J(t).  The  formalism  we 
develop  is  more  general,  however,  and  does  not  require  that  J(t)  be  Markov. 
Sytppose  that  B>(t)  is  governed  by  the  set  of  hazard  rates 
{v,.  i  a  •>  :  j,ke  J;  n,meW}.  Of  interest  in  this  paper  are  irreducible 
finite  Markov  chains  B(t)  for  which  N(t)  is  skip-free  in  both  directions 
so  that 


V(j,m),(k.n) 


0  if  |n-m|  >  1. 


(1.1) 


It  will  be  convenient  to  work  with  the  set  of  states  {(j,n)>  with 
common  row  index  n  as  an  entity,  and  to  introduce  the  corresponding 
notation  v+,  v°,  \f  to  designate  the  transition  rate  matrices  of  order 

*■7)  — n 

J+l  by 


^(j.iO.Oc.n)1 


(1.2) 


(b)  «  fv(j,n),k,n+l)^ 

(c)  «  ^V(j ,n) , (k,n-l)^  * 


We  will  also  work  with  the  matrix 

o  -  ♦ 

v  *  v  ♦  v  ♦  v 


(1.3) 


and  the  diagonal  matrix 


6 


v^n  -  diag  l  v(j>n)#(k>m)] 


(1.4) 


In  the  same  spirit  we  will  employ  the  transition  probability  matrices 


of  order  J+l 


2-(t)  ’  lp(j.»)(k,n)(t)1  ;  P(j.»),(k,n)(t)  *  P[B(t)  ■  (k.n)lB(O)  -  «,»)) 


(1.5) 


and  state  probability  row  vectors 


£„(t)  *  [P(j>n)(t))  ;  P(j>n)(t)  *  P[B(t)  «  (j,n)].  (1.6) 


The  ergodic  row  vectors  will  be  designated  by 


(1.7) 


where  e^  *  lim  p^  n^(t).  Laplace  transforms  will  often  be  used,  with 
the  notation  typical  of  subsequent  usage. 


I„(s)  *  L[p^(t)]  ■  /  £^(t)e’St  dt. 

We  will  also  employ  the  notion  of  a  matrix  p.d.f.  [5]. 


Def.  1.9  A  matrix  function  £(r)  *  [f__(r)l  is  a  matrix  p.d.f.  iff 

— — —  ran 


(1.8) 


a)  f  (x)  >  0 

UU1  — 


'm,n,x 


b)  l  J  f,k(x)dx  •  1  0<  j  <  J. 


k  -» 


Such  matrices  play  an  important  role  in  processes  defined  on  a  Markov 

chain  [7]  and  in  Markov  renewal  processes.  In  our  setting  f  (x)  ■  0  for 

ran 

x  <  0. 


i 


Examples 

To  motivate  the  analysis  that  follows  and  indicate  the  prevalence  of 
such  row- continuous  processes  some  examples  are  appropriate. 

A.  Contiguous  Processes 

Any  chain  B(t)  =  [J(t),N(t)J  on  8  ■  {(j,n)  |  0  <  j  <  J,  0  <  n  <  N} 
for  which  nj  (k  m)  *  0  if  |j-k|  >  1  or  |m-n|  >  1  will  be  called  con¬ 

tiguous.  The  row-continuous  chains  are  more  general  in  that  the  marginal 
row  process  need  not  be  skip-free. 

For  all  such  contiguous  processes  the  transition  rate  matrices  (1.2) 
are  tri-diagonal  (cf  Fig.  l.A).  All  such  processes  that  are  irreducible  are 
amenable  to  the  methods  we  describe. 

B.  Contiguous  Horizontal -Vertical  Processes 

A  subset  of  the  contiguous  processes  are  those  for  which  either  J(t) 
or  N(t)  can  change,  but  not  both  simultaneously.  If,  for  example,  J(t) 
and  N(t)  were  independent  truncated  birth-death  processes,  B(t)  would  be 
horizontal -vertical,  since  the  probability  of  simultaneous  change  would 
be  zero. 

Another  set  of  processes  ji(t)  *  [J (t) ,N(t) ]  has  J(t)  an  independent 
truncated  birth-death  process  and  N(t)  a  dependent  birth-death- like  process 
for  which  the  upward  and  downward  transition  rates  change  when  J(t)  changes. 
One  example  of  this  type  is  a  communications  link  carrying  both  voice  and 
data  [1].  We  note  that  for  such  processes  v*  and  v~  are  diagonal.  See 
Fig, l.B. 

C.  Tandem  Queues  with  Blocking 


A  contiguous  process  of  interest  is  the  tandem  two  station  series 
queue  where  each  station  has  finite  waiting  room,  arrivals  are  poisson  with 
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rate  A,  and  service  times  are  exponential  with  service  rates  n2-  For 
this  process  is  diagonal,  is  upper  diagonal,  and  is  lower 
diagonal.  (A  matrix  a  »  [a^]  is  upper  diagonal  if  a„  0  implies 
j  =  i+1.)  See  Fig.  l.C. 

The  methodology  treated  here  allows  each  station  to  consist  of  a 
finite  number  of  servers.  Various  types  of  blocking  and  feedback  could 
as  easily  be  analyzed,  although  the  example  of  §7  has  blocking  defined  by, 
the  first  queue  stops  serving  while  the  second  waiting  room  is  full. 

Extensive  numerical  analysis  and  further  discussion  of  this  tandem 
queue  model  can  be  found  in  §7. 

D.  Assembly  Line 

One  interesting  non-contiguous  process  is  an  assembly  line  with  two 
machines  in  sequence  and  finite  buffer  storage  between.  The  marginal 
process  J(t)  describes  fluctuating  working  and  non-working  states  of  the 
two  machines,  and  M2  which  are  governed  by  failure  rates  u^,  an(* 
repair  rates  A^,  A2  respectively.  The  two  machines  process  work  at  the 
same  rate  of  speed,  i.e.  have  a  common  hazard  rate  p  for  completing  their 
operation.  When  the  second  machine  is  down,  items  accumulate  in  the  buffer. 
When  the  first  machine  is  down,  the  flow  of  incoming  items  is  cut  off  and 
the  second  machine  goes  idle  after  the  buffer  is  emptied.  The  first  machine 
stops  when  the  buffer  is  full.  N(t)  is  the  number  of  half-finished  items 
in  the  buffer. 

For  this  model,  both  v*  and  v  "  are  diagonal  with  two  non-zero 

— n  »  n 

elements.  The  transition  rate  matrix  v  0  has  a  zero  diagonal,  and  (v°)T 

™  n  bH 

has  a  zero  diagonal.  See  Fig.  l.D. 
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(D)  Assembly  line  -  Non -contiguous 
(full  row  shown) 


Fig.  1.1  Possible  transitions  £rom  an  interior  point 
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§2.  Passage  Tine  Densities 

The  description  of  the  ergodic  and  dynamical  behavior  we  will  develop 
is  based  upon  the  passage  time  densities  for  the  bivariate  process.  To 
exploit  the  skip-free  feature  of  row- continuous  processes,  the  passage 
times  between  states  of  adjacent  rows  are  treated  as  building  blocks.  The 
focus  on  a  row  of  states  as  an  entity  then  gives  rise  to  a  matrix  probability 
density  function  describing  the  joint  distribution  between  the  time  of 

arrival  at  the  adjacent  row  and  the  state  reached,  as  a  function  of  state 

of  origin.  From  the  passage  time  p.d.f. 's  for  adjacent  rows  one  can  then 
describe  the  regeneration  times  for  the  states  of  a  row  and  hence  the 
ergodic  probabilities.  This  will  be  clearer  from  the  details  of  the 
development. 

To  discuss  the  passage  time  densities  it  will  be  useful  to  work  with 
the  lossy  process  J*(t)  for  row  n  on  the  set  {(j,n)>  for  which  other  states 
of  B  (i.e.,  other  rows)  are  absorbing.  Let 

P*(j,n),(k,n)(t)  =  P[J(t)  *  k’  N(t'}  “  0  <  *’  <  t  |  J(0)  =  j.  N(0)  =  n] 

and  as  for  (1.5)  let  p^^(t)  *  [p*^  ^  nj(t)].  Then  since 

ft  l$n(t)  a  -  fi^^Dn  +  we  rewrite  j£n(t)  in  the  form 

p*  (t)  »  exp{Q*t}  ,  (2.1) 

“hn  “n 

where  (P*  is  the  £  matrix  for  the  lossy  process,  i.e., 

£  -  -  V_  ♦  v°  .  (2.2) 

™T1  *Un 

Note  that  jg*  (t)  *s  strictly  substochastic.  From  (2.1)  one  then  has  for 
it*  (s)  «  L[p*  (t)] 

■Tin  “Tin 


,  s  >  0. 


(2.3) 


tt*  (s)  *  /  e'st  p*  (t)dt  «  [si  -  Q.*]'1 

”nn  o  *roi  — 


The  chain  B(t)  is  uniformizable  [6  ]  in  the  sense  that  there  exists 

a  v  such  that  “  >  v  >  max  l  v,.  .  .  when  J  and  K  are  finite, 

(j ,m)  (k,n)  CJ , (k,n) 

that  is  B  is  finite.  If  we  let 


♦ 

a 

-41 


h  i 


i  jS 


(2.3a) 


we  have  from  (2.3) 


tt*  (s) 
-nn 


l(s*v)£  -  va^]'1  . 


(2.4) 


To  proceed  further,  we  require  the  passage  time  densities  and  some 
associated  notation.  Let  s*.  (t)  be  the  joint  probability  density  that 

the  set  {(j,n*T)  :  jeJ)  is  visited  for  the  first  time  at  (k,n+l)  and  that 


the  time  of  visit  is  T,  so  that  J  /  s+, (t)dx  =  1.  Similarly,  let  s~... 

k  0  n»JK 

be  defined  with  respect  to  visiting  the  set  {(j,n-l)  :  jeJ}  for  n  ^  1.  The 

irreducibility  of  B(t)  implies  that  (cf  def.  1.9)  s+(t)  =  [s+  . ,  (t) }  and 

— n  n;jK 

fi^(t)  *  [s”.^(t))  are  matrix  p.d.f.  •  s.  Correspondingly, 

00  00 

*  /  j^ftidt  and  •  j  s^(t)dt  are  stochastic  matrices.  In  terms  of 


(t) 


the  Laplace  transforms,  c£(s)  *  L [s^ Ct) 3  and  cT(s)  *  L[s^(t)],  and  moment 

matrices  u*  *  f  t  s+(t)dt,  p"  *  /*  t  s*(t)dt  one  has 

-n  Q  -n  -n  q  -n 
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The  skip-free  feature  of  the  row- continuous  processes  provides  a 
recursion  between  the  upward  passage  time  p.d.f. 's  $*,  and  between  the 
downward  passage  time  p.d.f.'s  s^,  that  is  a  direct  matTix  counterpart  to 
that  present  for  birth-death  processes  where  J  *  {0}.  The  recursion  is 
based  on  an  argument  similar  to  that  employed  there  [6,4].  For  N(t)  to  go 
from  n  to  n+1,  it  must  either  do  so  directly  after  motion  on  row  n,  or 
there  must  be  a  first  downward  transition  to  row  n-1,  a  first  subsequent 
return  to  row  n,  and  then  a  first  subsequent  arrival  at  row  n+1.  A 
probabilistic  argument  based  on  this  then  gives  our  basic  recursive  equation 


t!(t) 


— n 


P*  +  C  \  1 

inn  — n  «nn  «-«n- 1  -41 


n  >  1 


(2.6) 


where  the  asterisk  denotes  convolution  in  time.  That  is,  £(t)*k(t)  * 

/*  a(t-T)b(x)dT.  Consequently,  from  (2.6)  and  (2.3a), 

n  “ 

o+(s)  *  vii*  (s)a+  ♦  vn*  (s)a"o+  (s)o+(s),  n  >  1  (2.7) 

—  n  «n  — n  -«n  1  —  n  ~ 

If  we  solve  (2.7)  for  o^(s)  and  use  (2.4)  we  obtain 


ii  -  *  *  »* 


-1 


(S)}  ]  o  (s) 


(rSJr:)  a* 

S+V  -n 


(2.8) 


In  place  of  (2.8)  one  has  for  n*0,  from  SpW  *  £jJ0(t)v0  and  (2.4) 


II  ‘ 


(-^— )  a* . 
vs+v‘,-o 


(2.9) 


Equations  (2.8)  and  (2.9)  can  be  used  to  generate  <£(s)  recursively.  Similar 
equations  generate  o"(s)  recursively  from  £^(s)  (cf  4.2).  By  letting  s«0 


we  get: 


{a0  ♦  a"  s*  }]s+ 
-u  -n  =n-l  —n 


o.  +  + 

a  Is  *  a 

*-0  “O  -o 


n  >  0 


(2.10) 


(2.11) 


two  recursion  relations  for  calculating  s  .  The  equivalent  relations  for 

— n 

s~  appear  in  (4.4). 

-n 

Differentiation  of  (2.8)  and  (2.9)  with  respect  to  s  at  s*0  leads  to 
the  following  relations  for  the  mean  first  passage  times: 


„  O  -  +  ,  ♦ 

[I  -  a  -  a  s  ,]y 
=■  — n  —n- 1  — n 


[I  -  a°]y+ 

“  “O  “O 


(2.12) 


(2.13) 


We  note  that  s*  is  stochastic  and  a0  ♦  a"  ♦  a*  is  stochastic.  Hence,  if 
— n  —n  — n  —  n 

is  not  zero  (guaranteed  by  irreducibility)  s^_j  is  strictly 

substochastic,  its  spectral  radius  is  less  than  one,  and  -  a°  -  a^  gt  .] 
is  invertible.  One  may  therefore  obtain  s^  from  (2.10)  and  from  (2.23). 
We  also  have  that  0^  <  o^(s)  *  for  s>®  real»  hence  (with  v/s+v  <  1)  we 
find  that  [^  -  (~j)  (a®  +  a^  *s  Invertible  for  s  ^  0  real. 

Therefore  (2.8)  may  be  used  to  obtain  ^J*(s) . 

Of  frequent  interest  in  applications  is  the  matrix  passage  time  p.d.f. 
s^^T)  describing  the  joint  distribution  of  the  time  at  which  row  n  is 
first  reached  and  the  state  reached  given  a  start  at  row  0.  Specifically 
s^^t)  ■  lson.ij(T)l  where  son.ij(T)  *  the  joint  probability  density  that 
the  set  {(j,n):jej)  is  visited  for  the  first  time  at  (j,n)  and  the  time 
of  first  arrival  is  t,  given  start  at  (0,i).  A  simple  probabilistic 
argument  shows  that  j^Ct)  -  iJ(t)*iJ(t)*...*jJ.1(t).  Correspondingly, 
2on(s)  *  £(»)£!<•).• -Cl(s>* 
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3.  Ergodic  Probabilities 

We  can  now  use  the  passage  time  densities  of  Section  2  to  find  the 
ergodic  probabilities,  as  outlined  in  the  introduction  to  that  section. 

The  basic  probabilistic  argument  for  finding  the  transition  probabilities 
from  the  passage  times  goes  as  follows.  If  one  is  in  row  n  at  time  0, 
to  be  in  row  n  at  time  t,  either 

a)  one  never  left  row  n 

b)  one  went  to  row  n+1  at  some  first  time,  returned  to  row  n  for 
the  first  time  subsequently,  and  was  in  row  n  at  time  t,  possibly  after 
further  wanderings, 

c)  one  went  to  row  n-1  at  some  first  time,  etc.  as  in  b) . 
Consequently,  one  has,  for  the  cases  0  <  n  <  N,  n«0,  and  n*N  respectively. 


EooCt)  ’ 

■  3k(t)  * 


From  (3.1) 


*(s) 

*rm 


it*  (s) 
-nn 


vit*  (s)a 

— nn  — 


A.i(s) 


"  (*> 
Tin 


(3.1) 


(3.1a) 


(3.1b) 


(3.2) 


•  ♦ 


♦  vjt*  (s )  a  a  (s)  jr  (s) 
Tin  ti  Ti  - 1  Tin 


,  0  <  n  <  N 


so  that 


+  . 


-  ♦ 


£„(>) 

Tin 


(3.3) 


One  then  finds,  from  (2.4),  that 


I  (s) 
nn 


(— — )[I  -  (— — )  {a0  *  a*  o'  (s)  +  a”  o+  (s)}l 
vs-*Vl-  Vs+v'  Si  -h  “h+1  “h  ^-1J 
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(3.4) 


To  show  that  the  matrix  expressed  within  brackets  is  in  fact  invertible 


we  let 


4<s, 


(— )  {a0  ♦  a+  o'  (s) 
s+v'  -h  -h  -h+l 


♦  a"  a*  (s)} 

-n  m-i 


n 


1 , . . . ,N-1 .  (3.5) 


o  ♦ 

Note  that  is  stochastic,  and  further 

L_1[8  (s) ]  «  (ve'VtI)*(a°  6(t)  ^  a+  s'  (t)  +  a'  s*(t)),  (3.6) 

=ti  ■“  *m  -n  "n+i  -^i  -n-1 


.  therefore  each  row  of  0  (s)  is  a  convex  combination  of  rows  of  matrix 

-•*n 

p.d.f.'s.  Hence,  B  (s)  is  itself  the  Laplace  transform  of  some  matrix 

■*n 

p.d.f.,  say  b  (t) .  Hence  fl-B  (s)]  is  invertible  for  s>0,  real.  With  this 
wi  *■  “n 

identification,  (3.4)  becomes 


2nn(s)  *  tl-Sn^r1  n-1.... , N-1.  (3.7) 

In  fact,  with  the  definition 

{£  *  £  £i(s)}  (3-8a) 

and 

6^(5)  *  {a^  ♦  a^  oN_i(s)}  (3.8b) 


we  immediately  get  (3.7)  for  all  n eW. 

T  T 

To  find  e  we  must  evaluate  lim  s  it  (s)  ■  1  e  .  In  [5]  it  is  shown  that 

“n  s-*Q+  “™  “  ^ 
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lin  s  ^(s)  •  i  •  -j - - 

«-m+  v  II 


*bn  iibnl 


,  T 
•  1  e. 

—  ^bn 


(3.9) 


where  e^  is  the  stationary  left  eigenvector  for  j^CO)  and  ubn  is  the  mean 


of  b^t),  that  is 


'  *bn 


Ubn  ‘  /„  1  V'>dt' 


(3.10) 


(3.11) 


We  can  easily  find  jg^O)  by  setting  s»0  in  (3.5),  (3.8a)  and  (3.8b),  while 
W^n  is  accessible  by  differentiation  via 

Ubn  *-&«»  '  *&&.!  ♦&&-!  "-1 . N-1 


i!bo  *  +  So  Si  ;  Sbn  *  v  +  Sn  i±N-l  .  (3.] 


As  we  will  see  in  Section  5,  row  balance 


^n+1  — n+1  — 


®T  &  I  »  0  <  n  <  N-l 


(3.14) 


is  always  present.  This  enables  one  to  evaluate 


"bn  ~  V  ^bn  %n  — 


(3.15) 


recursively  without  computing  u^n,  and  thus  and  jj^,  in  the  following 

T  IT 

manner.  From  (3.9)  one  has  e  ■  = —  e  so  that 

-n  ">bn  -  bn 


eT  ,  a'  1  -  - L- 

“n+1  -n*l  ~  mb>n4l 


ej  1  *'  ,  1  ;  ej  a+  1  ■  — —  eJ  a*  1 
-o,n^l  -4i +1  —  -n  -n  “  mb  n  _*>n  ■*  ~ 

(3.16) 
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Hence,  from  (3.15)  one  sees  that 


Vn+1 


*n+l  “ 

t — - 

^b,n  ah  — 


0  <  n  <  N-l  . 


(3.17) 


One  may  start  with  an  arbitrary  positive  M^0  and  then  normalize  at  the 
end.  More  specifically. 


T 


and 

K 


n^o  >  0,  arbitrary 


[  I 

ncN 


(3.18) 


(3.19) 


-19- 

4.  Summary  of  Computation  Procedure 

A  tabulation  of  the  key  results  obtained  above  is  given  to  provide 
an  overview  of  the  formalism,  and  ease  of  access  to  the  formulae  needed 
for  implementation. 


«  [I- 


*♦ 


(-54  a0]'1  (-24  a 

Vs*v  — 0J  s+v  — 0 


(4.1) 


£(S) 

— n 


[I  -  (— — • )  (a°  ♦  a*  a*  (s)  J1  ^  (-  — )  a* 

L-  Vs+v'  l— n  -n  — n-lv  J 11  -n 


Oj,(s)  c  l£  ■  (s^v) 


(4.2) 


o-(s) 

-41 


fl  -  (a0  ♦  A  a  fs)}]-1  (-^-)  a~ 

lmm  VS+V  -41  — n  -41  +  1  VS«V  -41 


♦ 

s 

-4) 


u  -  a°r‘  .* 

Vmm  -4>J  4) 


(4.3) 


♦ 

s 

-41 


[I  -{a°  ♦  a”  s*  .J]"1  a* 

—4i  -4i  —n-l  — 


-41 


ii  -  iS)'1  & 


(4.4) 


s 

-41 


II  -{£  *  &  i.i*!'1  '' 


In 


♦  1  n  0. -1  ♦ 

£o  v  ~  — oJ  -o 


(4.5) 


U+  ■  —  [I  -{a°  ♦  a*  s*  ,  }]”*  [I  ♦  va”  y+  .]  s+ 
in  v  —  in  -n-1  1  l—  ji  £n-lJ  ^ 


Ha 


v 


&  *  £  U ‘<s!!  ♦  &  i.i^"1 


(4.6) 


6  (S)  *  (— -)  [a  +  a  a*  (s)  ♦  a*  a  (s)] 

J  '“S+v  — n  -a  — n+1  *hti-1 

_  o  ♦  -  ♦ 

y0'  ■ 

v>  ■  fjjv>  il-  y*>rl 

4nil0)  ’  4n 

v  ■  \ ..I'iiVii yii>  ■  v>0 
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15.  Row  Balance  and  Row  Generation 


Proof : 

The  forward  equations  are: 

3t  *  Cct)  2$;  *  Cl  Cl  *  Cl(t)  Cl  CS-3) 

1  <n  <N-1 
and 

j£(t)  *  -VE^Ct)  ♦  «j  (5.3a) 

When  we  let  t-*«,  so  that  j£(t)  ♦  ej,  ^-^(t)  *  0T  the  result  is  immediate.  q 
Theorem  5.4 

If  a"  is  invertible  for  all  1  <  n  <  N  then 
■n  —  ~ 
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eT  =  {eT  [I  -  a°]  -  eT  a+  }(a\)_1  1  <  n  <  N-l 

-n+1  l-n  —  — nJ  -n-l  — n+i  -n+1  —  — 


(5.5) 


and 


4  -  4  £  *  £  tt)’1 


(5.5a) 


Proof :  trivial  by  Lenina  5.1.  g 

We  see  that  when  the  ergodic  distribution  is  desired,  and  this 

T 

invertibility  present,  the  e^  are  available  recursively  via  (5.4).  This 

T 

corollary  implies  that  one  need  calculate  only  one  eigenvector  (e^Q)  in 

T 

order  to  obtain  the  other  e  recursively. 


A  similar  result  holds  when  the  a*  for  0  <  n  <  N  are  invertible.  The 

— 

T 

recursion  then  begins  with  e^. 

Theorem  5.6 

If  a  is  invertible  for  all  n>0  then 
—  n 


T  ,  T  fT  o. 

e  ,  =  {e  [I  -  a  1 

-n-l  -n  —  — ir 


^n+1  in+1 


>  ‘i-P 


-1 


1, 


N-l 


and 

T 

^N-l 


(5.7) 


(5.7a) 


Proof : 

(5.7)  is  just  (5.2)  rearranged.  The  forward  equation  (5.3a)  has  its 
counterpart  on  the  top  row 

h  eJ(0  “  -  VeJ(0  ♦  VeI  2S1  +  VE^-ift)  fi.i  •  (5.8) 

When  t-**°  (5.8)  becomes  (5.7a) g 


Even  when  neither  set  is  invertible,  set  balance  can  be  used  to 
T 

normalize  the  «j)n*  as  described  at  the  end  of  section  3.  We  are  now  ready 
to  prove  our  basic  result. 


Theorem  5.9  (Row  balance) 

e  .  a  1  =  e  a  1 
-n+1  _n+i  —  — n  _n  — 

for  0  £  n  <  N-l 


(5.10) 


Proof : 

Recall  that  *  1_  -  (a^  +  1_  and  a^  1_  *  1_  -  a^  1_.  Combining 

this  with  Lemma  5.1  we  get 


T  ♦  .  T  -  , 

e  a„  1  -  e  a_.,  1 
-n  —  -n  -wi+1  — 


T  ♦  T 

e  .  a,  i  1  -  e*  a,  1 
-n-l  -n-l  —  -n  -wi  — 


(5.11) 


and 

eT  at  1  *  e?  a  7  1  (5.11a) 

— O  —  — 1  —1  — 


The  proposition  now  follows  inductively,  q 


i 

♦ 

i 
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§  6*  The  Matrix  Polynomial  Representation  of  2OJJ(s) 

Important  information  on  the  dynamical  behavior  of  a  Markov  chain 

is  contained  in  the  spectral  structure  of  its  first  passage  time  densities. 

Knowledge  of  this  structure,  and  that  of  the  related  relaxation  time,  is 

essential  when  one  wishes  to  use  the  ergodic  distribution  [6] . 

Towards  this  end  we  introduce  a  representation  of  a  (s)  of  the 

— mn 

form  Q^sJC^is),  where  Q^s)  is  a  matrix  polynomial. 

The  set  of  poles  of  o^s)  correspond  to  spectral  lines. 

Our  approach  closely  parallels  the  work  done  on  one-dimensional 
birth-death  processes  in  [2,3]  and  [6],  This  representation  will  be 
used  to  discuss  the  structure  of  s  (t)  and  simple  related  results.  We 
also  indicate  how  this  representation  may.be  used,  in  principle,  to  obtain 
the  relaxation  times. 

Thoerem  6.1 

If  a+  is  non-singular  for  all  n,  define 

—n 

2n.l(s)  ■  (>;)■'  HlTU-iDlW-itiW]  «>-2> 

with 

2c>(s)  *  I  ;  2l(s)  =  IC~)  I  -•£)•  (6.2a) 

Then 

the  matrix  polynomial  ^(s)  is  invertible  Vs  >  0  (6.3) 

*  2n<s>2;Vs’ 


(6.4) 
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Proof:  (By  induction) 


Clearly  (6.3)  is  true  for  n*0.  If  we  rewrite  O-j(s)  as  Qj(s)  « 

(~r)  (a*)  1  [I  -  (r~)  a0]  ;  a0  is  strictly  substochastic  hence 
Qj(s)  is  invertible  for  all  s  >  0.  Thus  (6.4)  makes  sense  and  is  clear¬ 
ly  true  for  n*0.  Now  assume  Q  (s)  is  non-singular  for  0  <  n  <  M  and 

•*1  —  — 

2>>  -  (s)  for  0  <  n  <  M-l.  Then 


W>  '  i  -  i£p V5)  -  Vi(s» 


(6.5) 


■  II-  ti&V  -  Si4i(5),12«w 


By  the  same  argument  as  that  for  (2.8)  et  al.,  we  get  that  ^ 

is  invertible  V  s  >  0.  Postmultiplying  (6.2)  by  Q^fs)  and  inverting 
—  -^n 


one  has  *  V^Cs),  fro.  (2.8) 


We  note  immediately  that  the  recursion  relation  (6.2)  implies  that 
V  s))^  is  a  polynomial  in  s  of  degree  n.  The  decomposition  (6.4)  will 
be  seen  to  be  useful  in  a  variety  of  ways. 

The  g^fs)  arrays  allow  one  to  evaluate  first  passage  times  upwards 
and  downwards  over  a  number  of  rows  rather  easily.  The  following  theorems 
illustrates  this: 

Theorem  6.6 

In  the  context  of  Theorem  6.1, 

.-1 


a)  2on(s)  •  (s)  for  n  >  1 


•  (<)■'  < i  v0’  ♦  a  -  iX<°>  -  > 


with  9^(0)  •  Oand^(O)  -  £  (a^) 


♦.-1 


b)  Q*+1(0) 


c>  iion  *  & 


Proof : 


a)  o  (s) 

~*on 


V  £m  ■  2,<s>  £1c*>s.(*»s;1(*> 

o*0  m*l 


,-l 


2n  (s)  by  Theorem  6.1 


b)  From  (6.2)  we  get 


&.!<*>  *  (&)'1  <j2„C5>  ♦  aX") 


(6.7) 


At  s*0  this  becomes  (b) .  The  cases  n*0,l  are  trivial. 


c)  From  (a)  2n^s^on^  *  JL»  hence 


&(s)Oon(s)  ♦  2n(s)oin(s)  *  0. 


(6.8) 


Therefore , 

not,  ■  -  ai„(°)  *  a'toJUltojai'to)  0  (6.9) 

Theorem  6.6  provides  a  computationally  easy  way  to  calculate  the 
mean  upwards  passage  times  when  the  a*  are  invertible.  A  dual  argument 
eases  the  computation  of  downwards  passage  times  when  the  aj^  are  invertible. 
The  explicit  formulation  for  arbitrary  upwards  passage  times  follows. 
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Theorem  6.10 


>-l 


a)  2inn^  *  (s)  for  0  <  m  <  n  <  N 


b)  J4m  =  V0)S1(0)^n(0)i1(0)  '  Si<°>Sil(°)  for  0  <  m  <  n  <  N 


Proof : 


a)  Since  a  (s)o  (s) 
— om  -mn 


a  (s)  the  result  follows  from  (6.6a) 
-on  v 


b)  We  use  (a)  to  obtain  P  (s)^(s)  *  differentiating  at  s=0, 

one  finds  that 


o*  (0)0  (0)  ♦  o  (0)2* (0)  =  OHO). 

™Tnn  ™ ti  -Tim  — n 


(6.11) 


From  (a)  and  the  identity  =  -  0^(0)  statement  (b)  follows.^ 


-1 


We  note  that  fs)  x(s)  implies  that 


-1 


(6.12) 


Thus,  calculation  of  (^(0) ,  and  2^(0)  is  sufficient  to  give  s^  and 
directly  (i.e.  non-recursively,  and  efficiently). 

The  matrices  2^(0)  and  (^(0)  can  be  calculated  recursively  using 
(6.2)  at  s*0,  and  (6.6b).  Knowledge  of  the  matrices  allows  us  to 
evaluate  the  mean  ergodic  exit  time  and  mean  stationary  sojourn  time, 
two  useful  measures  of  the  dynamical  behavior  of  the  queue  [6] .  Both  are 
defined  with  respect  to  a  partition  of  the  state  space  into  two  disjoint 
connected  sets,  called  the  good  set  and  the  bad  set. 


The  ergodic  exit  time  is  the  time  required  to  leave  the  good  set 
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given  that  the  system  has  settled  down,  i.e.  is  at  ergodicity  (and,  of 
course,  that  it  is  in  the  good  set).  Thus,  the  mean  ergodic  exit  time. 


Gm  *  {(j,n)  :  n  <  m}  ,  Bm 


is  given  by 


*  Ergodic  exit  time  from  G 


{(j,n)  :  n  >  m) 


(6.13) 


(6.14) 


^Vim  e^  1 

—  j  n<n  -n  — 


since  entry  into  is  at  row  m.  Using  Theorem  6.10,  we  get 

l  J{Q  (0)Q'1(0)Q'(0)0‘1(0)  -  Q'WQ'^O)}  1 
n<m  ^  n  -n  -n 

UEm  =  J  ~  :  * 

“el 
n<m  -n  — 


(6.15) 


The  stationary  sojourn  time  is  the  time  required  to  first  leave  the 
good  set  given  that  the  system  was  stationary  and  a  transition  into  the 
good  set  from  the  bad  set  just  occurred.  We  have 


=  Stationary  sojourn  time  on  G 


/  T  - 
(e  a 

-m  -4D 


Jk.  1  i}/ 


(1  -  , 
e  a  1 
—m  — 


(6.16) 


This  can  be  rewritten  in  a  simpler  form  using  a  well  known  result  [6]  as 


^Vm  * 

r  IB  tn 


(6.17) 
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§7.  A  Tandem  Queue  Example 

A  tandem  queue  is  a  queueing  system  with  distinct  service  facilities 
connected  in  series,  i.e.  the  customer  output  stream  of  the  first  facility 
is  the  input  stream  of  the  second  [9].  To  illustrate  the  algorithmic 
procedure  we  have  developed,  one  such  tandem  queue  will  be  analyzed. 

The  example  selected  has  been  chosen  primarily  for  didactic  reasons. 

A  more  complex  and  realistic  example  could  have  been  analyzed  with  little 
increase  in  machine  cost. 

The  tandem  queue  evaluated  is  a  two  server  series  queue  with  block¬ 
ing  and  finite  buffers  (see  il.E).  The  first  service  facility  has  8 
buffer  slots  and  4  servers.  The  second  facility  has  4  buffer  slots  and 
5  servers.  A  flow  diagram  is  given  in  Figure  7.1.  The  corresponding 
rate  matrices  are  given  in  Figure  7.2.  When  the  first  queue  is  full 
customers  balk  and  go  elsewhere.  When  the  second  queue  is  full,  the 
first  queue  stops  serving  until  a  space  at  the  second  facility  opens  up. 
The  model  can  easily  be  modified  to  allow  different  types  of  blocking, 
and  features  such  as  feedforward,  feedback,  etc. 

In  Figure  7.1  the  occupancy  level  at  the  first  facility  is  given  by 
the  coordinate  j  which  corresponds  to  the  number  of  people  in  service  or 
waiting  there.  The  occupancy  level  at  the  second  facility  is  given  by 
the  coordinate  n.  The  blocking  may  be  seen  in  the  absence  of  the 
transition  rates  associated  with  on  the  top  row. 

Figure  7.2  displays  the  three  iratrices  \£,  and  when  n«3. 

The  matrix  is  upper  diagonal  (X)  for  all  n.  The  only  transitions 
which  have  no  impact  on  the  2nd  queue  are  the  arrivals  (at  rate  A)  to 
the  first  queue,  which  increase  that  occupation  level  by  1.  The  matrix 
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is  diagonal,  with  diagonal  elements  u^mnininyk^}  where  k£  is  the 

number  of  servers  in  the  2nd  service  facility.  Finally,  v+  is  lower 

-n 

diagonal.  Increases  in  queue  length  at  the  2nd  service  facility  are 

caused  by  departures  from  the  first  facility. 

The  results  are  shown  in  Figs.  7.3  (a) , (b) , (c) ,  7.4  and  7.5  for  a 

traffic  intensity  of  0.4  at  the  first  facility  and  0.6  at  the  second. 

Figs.  7.3  (a),(b),(c)  describe  the  ergodic  probabilities  as  a  function 

of  U,n).  An  examination  of  these  figures  inraediately  reveals  two  features. 

First,  the  ergodic  probabilities  are  signficant  for  modest  values  of  j 

and  n,  arising  from  the  moderate  traffic  intensities.  The  small  ridge 

visible  along  the  upper  edge  of  Fig.  7.3  (a)  indicates  some  blocking. 

The  ergodic  probabilities  may  be  useless  when  system  parameters 

change  before  ergodicity  is  reached.  The  information  contained  in  the 

first  passage  times  may  then  be  helpful.  The  ergodic  distribution,  for 

example,  might  cause  concern  over  large  probabilities  of  saturation, 

blocking,  long  waiting  times  and  so  on.  If,  however,  we  look  (cf.  Fig. 

7.5)  at  the  mean  first  passage  times  from  idleness,  we  see  that  the 

mean  time  to  saturation  is  on  the  order  of  4  days  (when  our  time  scale 

is  in  hours).  If  this  queue  modeled  a  system  that  starts  anew  each  day, 

one  would  be  less  inclined  to  worry. 

The  tables  in  Figure  7.4  give,  for  differing  levels,  the 

mean  sojourn  time  [6],  i.e.,  the  mean  time  spent  above  a  certain  level 

after  the  level  is  reached.  This  gives  a  feeling,  for  example,  for 

the  persistence  time  in  a  congested  state.  The  table  7.5  presents 


mean  first  passage  times  one  step  upwards  (downwards).  The  (j,n)  entry 
in  the  table  is  the  mean  time  to  go  from  (j,n)  to  a  state  in  the  row 
n+1  (n-1).  Again,  the  effect  of  blocking  is  noticeable  in  the  dropoff 
(for  the  upward  table)  of  the  passage  time  with  increasing  values  of  j . 

Ergodic  distributions  for  the  same  system  with  different  parameters 


are  shown  in  Figure  7.6 


for  comparison  purposes. 
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7.6  Ergodic  probabilities  of  tandem  queues  with  different  parameters 
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§8.  Rank  Reduction  in  the  Row-Continuous  Model 

The  algorithms  presented  in  the  preceding  sections  exploit  the 

block  tridiagonal  structure  of  the  transition  rate  matrix  associated 

with  a  bivariate  row- continuous  markov  chain.  These  algorithms  are 

considerably  more  efficient  than  the  general  algorithms  for  computation 

on  a  markov  chain  both  in  computer  time  needed  to  attain  desired  results 

and  in  storage  needed  to  perform  these  computations.  To  illustrate  this 

fact  we  compare  the  resources  required  to  calculate  ergodic  probabilities 

for  the  tandem  queue  of  section  7. 

We  note  that  the  tandem  queue  is  not  time  reversible,  and  so  we 

could  not  have  used  detailed  balance  simplification.  We  also  note  that 

the  mean  passage  times,  themselves  of  some  importance,  are  a  byproduct  of 

the  row-continuous  ergodic  distribution  calculations.  Finally,  the 

tandem  queue  does  have  the  property  that  is  trivially  invertible, 

a  feature  which  reduces  computation  time. 

In  order  to  calculate  ergodic  probabilities  efficiently  we  must 

store  the  a£.  of  course,  as  well  as  the  s*  and  s^.  The  storage 

2 

requirements  are  then  0(J  N) .  For  a  contiguous  problem  (such  as  this) 

one  would  pick  the  row  direction  J(t)  to  be  that  with  the  fewest  states. 

A  naive  generalized  approach  would  store  the  full  transition  matrix, 

2  2 

thus  using  up  0(J  N  )  memory  locations.  For  example,  one  tandem  queue 
problem  with  J*8  and  N*80  was  evaluated  in  a  workspace  of  140,000  bytes. 
The  transition  rate  matrix  for  the  bivariate  problem  would  require  more 
than  4,000,000  bytes  alone. 

The  example  of  section  7  was  coded  in  APL,  a  notoriously  inefficient 
language.  Table  8.1  gives  some  times  for  the  tandem  queue  evaluation. 
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These  were  run  on  an  Amdahl  470  v/6  timesharing  with  20  users.  The 
timings  include  all  setup  and  all  calculations  required  to  find  the 
ergodic  probabilities.  For  comparision  purposes,  a  single  eigenvalue 
was  found  for  a  64  x  64  transition  matrix.  This  required  approximately 
70  CPU  seconds  using  an  optimized  successive  approximation.  No  effort 
was  made  to  optimize  the  programs  used  for  Table  8.1.  In  fact,  the 
matrices  were  regenerated  when  needed. 

A  complexity  argument  can  be  used  to  find  the  number  of  operations 

required  to  evaluate  the  ergodic  distribution.  This  would  severely 

understate  the  advantage  of  using  the  previous  algorithms,  because  a 

general  approach  would  require  repeated  paging  in  and  out  of  the  system 

parameters  for  a  problem  of  any  size.  Nevertheless,  the  run-time  was 

2 

determined  empirically  as  0(J  N).  We  expect  that  when  it  is  necessary  to 
calculate  rather  than  being  able  to  explicitly  determine  them 

(as  in  the  tandem  queue  where  is  diagonal)  the  calculation  time  would 
be  0(J3N),  still  a  greyt  improvement  over  the  general  case. 
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J _ N _ CPU  time  (seconds  on  an  IBM  370-158) 
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